//
// Created by me1on on 22-5-18.
//
#define SOLAR_BUILD_CUDA_REFINED_CYLINDER_RECEIVER
#include "CudaRefinedCylinderReceiver.cuh"

#include <algorithm/rayTracing/qmcrt_util/global_constant.h>
#include <algorithm/rayTracing/qmcrt_util/global_function.cuh>
#include <algorithm/rayTracing/qmcrt_util/vector_arithmetic.cuh>

#include <cmath>
#include <cstdio>

using namespace solar::cuda;

void CudaRefinedCylinderReceiver::CInit(float pixel_per_meter_for_receiver)
{
    pixel_length_ = 1.0f / float(pixel_per_meter_for_receiver);
    Cset_resolution(pixel_per_meter_for_receiver);
    Calloc_image();
    Cclean_image_content();
}

void CudaRefinedCylinderReceiver::Cset_resolution(float pixel_per_meter_for_receiver)
{
    // resolution_.x = int(ceil(2 * M_PI * (M_PI/2 + asinf(radius_/size_.x)) * size_.x *
    // float(pixel_per_meter_for_receiver)));
    //    resolution_.x = int(ceil(M_PI * M_PI * size_.x * float(pixel_per_meter_for_receiver)));
    //    resolution_.y = int(size_.y * float(pixel_per_meter_for_receiver));

    //    resolution_.x=panel_*tube_*tubeUnit_;
    resolution_.x = panel_ * (tube_ * tubeUnit_ + tube_ - 1);
    //    resolution_.x=18*300;
    resolution_.y = ceil(size_.y * float(pixel_per_meter_for_receiver));
    printf("-------refined cylinder resolution = (%d, %d)\n", resolution_.x, resolution_.y);
    //    printf("image %d\n",resolution_.x*resolution_.y);
}

float3 CudaRefinedCylinderReceiver::getFocusCenter(const float3& heliostat_position)
{
    if (isInsideCylinder(heliostat_position))
    {
        // If the heliostat position is inside cylinder, return the receiver position
        return pos_;
    }

    float3 dir = normalize(heliostat_position - pos_);
    float radius = size_.x / (length(make_float2(dir.x, dir.z)));

    float x = pos_.x + dir.x * radius;
    float z = pos_.z + dir.z * radius;
    //    printf("focus %f %f %f\n",x,pos_.y,z);
    return make_float3(x, pos_.y, z);
}

__device__ __host__ auto
CudaRefinedCylinderReceiver::GIntersect(const float3& orig, const float3& dir, float& t, float& u,
                                        int& panel_no, int& tube_no, int& angle_no,
                                        int& ifinterspace, int& interspace_no) -> bool
{
    // If the origin is inside the cylinder, it won't intersect with it
    if (isInsideCylinder(orig))
    {
        // printf("The origin is inside the cylinder!!!\n");
        return false;
    }

    float r = radius_;
    float panel_angle = 2 * MATH_PI / float(panel_);

    float R = size_.x;
    float2 Ro = make_float2(pos_.x - orig.x, pos_.z - orig.z);
    float tp = dot(Ro, normalize(make_float2(dir.x, dir.z)));
    float delta = dot(Ro, Ro) - tp * tp;

    // Return false if:
    //  1) The direction is different
    //  2) Has no intersection
    //        float R2 = (R+extend) *(R+extend) ;   // radius ^ 2
    float R2 = (R) * (R); // radius ^ 2
    if (tp < -Epsilon || delta > R2)
    {
        return false;
    }

    // remain and calculate the positive solution
    float t_plus = delta <= 0.0F ? (R) : sqrt(R2 - delta);
    t = (tp - t_plus) / length(make_float2(dir.x, dir.z));

    // calculate intersection point
    float3 intersect_pos_cylinder = t * dir + orig;
    float h_rece = size_.y + 0.5;
    u = (intersect_pos_cylinder.y - pos_.y) / h_rece + 0.5f;
    if (u < 0.0F || u > 1.0F)
    {
        return false;
    }
    // intersect_pos_1 is to intersect between big cylinder and reflected ray
    float3 IO = intersect_pos_cylinder - make_float3(0, intersect_pos_cylinder.y, 0);
    // IO is the vector from center of the cylinder to intersect pos 1
    // theta_IO is the angle from (1,0,0) East to IO
    //  float3 east= make_float3(1,0,0);
    float3 north = make_float3(0, 0, -1);
    float theta_IO = 0;
    if (IO.x > 0)
    {
        theta_IO = dot(IO, north) / length(IO);
        theta_IO = acos(theta_IO);
    }
    else
    {
        theta_IO = dot(IO, north) / length(IO);
        theta_IO = 2 * MATH_PI - acos(theta_IO);
    }
    panel_no = int(floor(theta_IO / panel_angle)) % panel_;
    float panel_up_theta = float(panel_no) * panel_angle;
    float panel_down_theta = float(panel_no + 1.0) * panel_angle;

    float up_x = sin(panel_up_theta) * R;
    float up_z = -cos(panel_up_theta) * R;
    float down_x = sin(panel_down_theta) * R;
    float down_z = -cos(panel_down_theta) * R;

    float Height = pos_.y;
    float half_H = size_.y / 2;

    float3 panel_up = make_float3(up_x, Height + half_H, up_z);
    float3 panel_down = make_float3(down_x, Height + half_H, down_z);

    float3 up_down = normalize(panel_down - panel_up);

    float3 mid = make_float3((up_x + down_x) / 2, panel_up.y, (up_z + down_z) / 2);
    float3 normal_to_cylinder = make_float3(mid.x, 0, mid.z);
    normal_to_cylinder = normalize(normal_to_cylinder);
    float3 panel_C = make_float3(panel_up.x, panel_up.y - size_.y, panel_up.z);
    // Down   UP
    //
    //        Panel_C(x,z same with up, but y lower than up)
    float t_ = 0;
    float u_ = 0;
    float v_ = 0;
    float tube_unit_angle = MATH_PI / tubeUnit_;

    bool isIntersectWithMargin = global_func::rayParallelogramIntersect(
        orig, dir, panel_up, panel_down, panel_C, t_, u_, v_);
    if (isIntersectWithMargin == false)
    {
        panel_no = (panel_no + 1) % panel_;
        panel_up_theta = float(panel_no) * panel_angle;
        panel_down_theta = float(panel_no + 1) * panel_angle;

        up_x = cos(panel_up_theta) * R;
        up_z = sin(panel_up_theta) * R;
        down_x = cos(panel_down_theta) * R;
        down_z = sin(panel_down_theta) * R;

        panel_up = make_float3(up_x, Height + half_H, up_z);
        panel_down = make_float3(down_x, Height + half_H, down_z);
        mid = make_float3((up_x + down_x) / 2, panel_up.y, (up_z + down_z) / 2);
        normal_to_cylinder = make_float3(mid.x, 0, mid.z);
        normal_to_cylinder = normalize(normal_to_cylinder);
        panel_C = make_float3(panel_up.x, panel_up.y - size_.y, panel_up.z);

        isIntersectWithMargin = global_func::rayParallelogramIntersect(
            orig, dir, panel_up, panel_down, panel_C, t_, u_, v_);
        if (!isIntersectWithMargin)
        {
            panel_no = (panel_no - 2) % panel_;
            panel_up_theta = float(panel_no) * panel_angle;
            panel_down_theta = float(panel_no + 1) * panel_angle;
            up_x = cos(panel_up_theta) * R;
            up_z = sin(panel_up_theta) * R;
            down_x = cos(panel_down_theta) * R;
            down_z = sin(panel_down_theta) * R;
            panel_up = make_float3(up_x, Height + half_H, up_z);
            panel_down = make_float3(down_x, Height + half_H, down_z);
            mid = make_float3((up_x + down_x) / 2, panel_up.y, (up_z + down_z) / 2);
            normal_to_cylinder = make_float3(mid.x, 0, mid.z);
            normal_to_cylinder = normalize(normal_to_cylinder);
            panel_C = make_float3(panel_up.x, panel_up.y - size_.y, panel_up.z);

            isIntersectWithMargin = global_func::rayParallelogramIntersect(
                orig, dir, panel_up, panel_down, panel_C, t_, u_, v_);

            if (!isIntersectWithMargin)
            {

                return false;
            }
        }
        // 和内平板是否有交点  为啥和外平板没有交点 和内平板有交点？
        // 回答：外平板在前 内平板在后 有可能光线会打在内平板上下左右边缘
    }

    float3 intersect_with_panel = orig + t_ * dir;
    t = t_;
    // 交在内侧面上。说明肯定和圆柱有交点。
    // 根据面板上交点位置，计算与该圆柱交点 以及t
    // 循环往左往右寻找上一个相交的圆柱 以及t
    // 直到与上一个圆柱不相交，停止 与当前圆柱相交，计算交点，返回
    // dir intersect with panel
    //  next to judge the no of little cylinder
    // litter cylinder raidus  and
    float D = length(panel_up - panel_down);
    float AO = u_ * D;

    float interspace = D - float(tube_) * 2 * r;
    // assmue there N tube and N-1 interspace
    interspace = interspace / (float(tube_) - 1);
    int inter_tube = floor(AO / (radius_ * 2 + interspace));
    float diss = inter_tube * (radius_ * 2 + interspace);
    int LR = 0; // 往前还是往后找 前=1 后=2  前就是往up找 后就是往down找

    if (AO - diss >= 2 * r)
    {
        // 交点在interspce处
        // 求左 求右交点 如果都没有 交点就在interspace上
        // 否则
        float3 tube_pos = panel_up + ((2 * r + interspace) * inter_tube + r) * up_down;
        float2 r0 =
            make_float2(tube_pos.x - intersect_with_panel.x, tube_pos.z - intersect_with_panel.z);
        float tp0 = dot(r0, normalize(make_float2(-dir.x, -dir.z)));
        float delta_0 = dot(r0, r0) - tp0 * tp0;
        float R_2 = r * r;
        if (tp0 < -Epsilon || delta_0 > R_2)
        {
            float3 tube_posb = panel_up + ((2 * r + interspace) * (inter_tube + 1) + r) * up_down;
            float2 r0b = make_float2(tube_posb.x - intersect_with_panel.x,
                                     tube_posb.z - intersect_with_panel.z);
            float tp0b = dot(r0b, normalize(make_float2(-dir.x, -dir.z)));
            float delta_0b = dot(r0b, r0b) - tp0b * tp0b;
            if (tp0b < -Epsilon || delta_0b > R_2)
            {
                // 交点在空隙处，和左右两个小圆柱都没有交点
                ifinterspace = 1;
                interspace_no = inter_tube;
                u = v_;
                //                    printf("tp0b %f\n",tp0b);
                return true;
            }
            LR = 2;
        }
        else
        {
            LR = 1;
        }
    }
    float3 tube_pos_1 = panel_up + (diss + r) * up_down;

    float R_2 = r * r;
    // 肯定有交点
    float ttt_1 = radius_;
    // 如果在最前头，往后找
    // 如果在最后头，往前找
    if (inter_tube == 0)
    {
        LR = 2;
    }
    if (inter_tube == tube_ - 1)
    {
        LR = 1;
    }

    if (LR == 0)
    {
        // 还不知道往左还是往右找的，往左求，往右求
        // 前 inter_tube-1
        // 后 inter_tube+1
        // 先前
        int front_tube = inter_tube - 1;
        float ftdiss = front_tube * (radius_ * 2 + interspace);
        float3 fttube_pos_1 = panel_up + (ftdiss + r) * up_down;
        float2 ftr0_1 = make_float2(fttube_pos_1.x - intersect_with_panel.x,
                                    fttube_pos_1.z - intersect_with_panel.z);
        float fttp0_1 = dot(ftr0_1, normalize(make_float2(-dir.x, -dir.z)));
        float ftdelta_0 = dot(ftr0_1, ftr0_1) - fttp0_1 * fttp0_1;
        if (fttp0_1 < -Epsilon || ftdelta_0 > R_2)
        {
            int back_tube = inter_tube + 1;
            float btdiss = back_tube * (radius_ * 2 + interspace);
            float3 bttube_pos_1 = panel_up + (btdiss + r) * up_down;
            float2 btr0_1 = make_float2(bttube_pos_1.x - intersect_with_panel.x,
                                        bttube_pos_1.z - intersect_with_panel.z);
            float bttp0_1 = dot(btr0_1, normalize(make_float2(-dir.x, -dir.z)));
            float btdelta_0 = dot(btr0_1, btr0_1) - bttp0_1 * bttp0_1;
            if (bttp0_1 < -Epsilon || btdelta_0 > R_2)
            {
                // 如果这里确定了与左右都不交，可以返回这里
                //                    return false;
                float3 intersect_pos_tube = intersect_with_panel - ttt_1 * dir;

                u = (intersect_pos_tube.y - pos_.y) / size_.y + 0.5f;
                if (u < 0.0f || u > 1.0f)
                {
                    return false;
                }
                float3 intersect_to_centertube = intersect_pos_tube - tube_pos_1;
                intersect_to_centertube.y = 0;
                float intersect_thetatube =
                    acos(dot(intersect_to_centertube, -up_down) /
                         (length(up_down) * length(intersect_to_centertube)));
                if (intersect_thetatube > MATH_PI)
                {
                    return false;
                }
                tube_no = inter_tube;
                angle_no = floor(intersect_thetatube / tube_unit_angle);
                //                        angle_no=int(tubeUnit_/2)-1;
                t = t - ttt_1;
                return true;
            }

            LR = 2;
        }
        else
        {
            LR = 1;
        }
    }

    int flag = 1;

    int nowtube = inter_tube;
    float3 tubeposs = tube_pos_1;
    float tend = ttt_1;
    while (true)
    {
        int inter_tubeLR;
        if (LR == 1)
        {
            inter_tubeLR = nowtube - 1;
        }
        else
        {
            inter_tubeLR = nowtube + 1;
        }
        float dissLR = inter_tubeLR * (radius_ * 2 + interspace);
        float3 tube_pos_LR = panel_up + (dissLR + r) * up_down;
        float2 r0_1LR = make_float2(tube_pos_LR.x - intersect_with_panel.x,
                                    tube_pos_LR.z - intersect_with_panel.z);
        float tp0_1LR = dot(r0_1LR, normalize(make_float2(-dir.x, -dir.z)));
        float delta_0LR = dot(r0_1LR, r0_1LR) - tp0_1LR * tp0_1LR;
        if (tp0_1LR < -Epsilon || delta_0LR > R_2)
        {
            // 这里返回的是上一步的
            float3 intersect_pos_LR = intersect_with_panel - tend * dir;
            u = (intersect_pos_LR.y - pos_.y) / size_.y + 0.5f;
            if (u < 0.0f || u > 1.0f)
            {
                return false;
            }
            float3 intersect_to_centertube = intersect_pos_LR - tubeposs;
            intersect_to_centertube.y = 0;
            float intersect_thetatube = acos(dot(intersect_to_centertube, -up_down) /
                                             (length(up_down) * length(intersect_to_centertube)));
            if (intersect_thetatube > MATH_PI)
            {
                return false;
            }
            tube_no = nowtube;
            angle_no = floor(intersect_thetatube / tube_unit_angle);
            //                    angle_no=int(tubeUnit_/2)-1;

            t = t - tend;
            return true;
        }
        else
        {
            // 如果当前的到头了，返回当前的，否则继续循环

            float t_plus_LR = delta_0LR <= 0.0f ? r : sqrt(R_2 - delta_0LR);
            tubeposs = tube_pos_LR;
            tend = (tp0_1LR + t_plus_LR) / length(make_float2(dir.x, dir.z));
            nowtube = inter_tubeLR;

            if ((inter_tubeLR == 0 && LR == 1) || (inter_tubeLR == tube_ - 1 && LR == 2))
            {
                float3 intersect_pos_LR = intersect_with_panel - tend * dir;
                u = (intersect_pos_LR.y - pos_.y) / size_.y + 0.5f;
                if (u < 0.0f || u > 1.0f)
                {
                    return false;
                }
                float3 intersect_to_centertube = intersect_pos_LR - tubeposs;
                intersect_to_centertube.y = 0;
                float intersect_thetatube =
                    acos(dot(intersect_to_centertube, -up_down) /
                         (length(up_down) * length(intersect_to_centertube)));
                if (intersect_thetatube > MATH_PI)
                {
                    return false;
                }

                tube_no = nowtube;
                angle_no = std::floor(intersect_thetatube / tube_unit_angle);
                //                        angle_no=int(tubeUnit_/2)-1;

                t = t - tend;
                return true;
            }
        }
        flag++;
    }
}

__device__ __host__ auto CudaRefinedCylinderReceiver::isInsideCylinder(const float3& orig) -> bool
{
    float2 l = make_float2(orig.x - pos_.x, orig.z - pos_.z);
    return dot(l, l) <= size_.x * size_.x;
}
